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Optical homodyne tomography is discussed in the context of classical image pro- 
cessing. Analogies between these two fields are traced and used to formulate an 
iterative numerical algorithm for reconstructing the Wigner function from homo- 
dyne statistics. 
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1 Introduction 



Several years after the first demonstration [1], optical homodyne tomography has be- 
come a well established tool for measuring quantum statistical properties of optical 
radiation. What is particularly fascinating, this technique provides practical means 
to visualise the measured quantum state in the form of the Wigner function. This 
success is a result of combining a complete quantum mechanical measurement of field 
quadratures with a filtered back-projection algorithm used in medical imaging. 

The purpose of this contribution is to trace some other analogies between quantum 
state reconstruction and classical image processing, with the motivation to develop novel 
numerical methods for quantum tomography. Our interest will be focused on imper- 
fect detection [2, 3], which has deleterious effects on quantum interference phenomena 
exhibited by non-classical states [4]. As we will discuss in Sec. 2, such effects can be 
related in the phase space representation to image blurring. Restoration of blurred im- 
ages is a well known problem in classical imaging, and a number of methods has been 
developed for this purpose. Specifically, we shall briefly describe in Sec. 3 the Richard- 
son algorithm [5] (known also in statistics as the expectation-maximization algorithm 
[6]), which provides an effective iterative procedure to perform image deblurring. 
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An interesting question is, whether classical deblurring methods can be transferred 
to quantum tomography. We discuss this, in the case of the Richardson algorithm, in 
Sec. 4. The answer is not straightforward: the Richardson algorithm assumes positive 
definiteness of the original undegraded image, and this condition is not satisfied by 
the quantum mechanical Wigner function. We will show that this difficulty can be 
overcome by expressing the Wigner function in terms of the phase space displaced 
photon distribution. This yields an iterative algorithm for reconstructing the Wigner 
function, which incorporates compensation for detection losses in a numerically stable 
way [7]. 



2 Imperfect detection and image blurring 

Homodyne detection is a realization of the quantum mechanical measurement of field 
quadratures only in the idealized loss-free limit. In practice, a fraction of the signal field 
is always lost due to the mode-mismatch and the non-unit efficiency of photodiodes. 
The homodyne statistics collected using a realistic setup is described by the distribution 
[8]: 

h(x; 9) = r dx' 1 cxp f J x - Vvx') 2 \ (1) 
J-oo v/7r(l-r?) V 1 ~V J 

where 9 is the local oscillator phase, rj characterizes the overall detector efficiency, and 
(x'g\p\x' s ) denote diagonal elements of the density matrix p in the quadrature basis. 

Application of the back-projection transformation to h(x;9) yields, instead of the 
Wigner function, a generalized phase space quasidistribution function P(q,p; s) with the 
ordering parameter s = — (1 — 77) / 77. This function can be expressed as a convolution of 
the Wigner function W(q,p) with a gaussian function: 

P(q,p;-\s\) = Jdq'dp' ^exp(-±l( q - q r + (p-p r ) 2 ^W( q ',p>). (2) 

Thus what we reconstruct from imperfect homodyne data, is a blurred version of the 
Wigner function. The question is, whether we can get rid of this blurring in numerical 
processing of experimental data. 

A similar problem appears in the following classical context: suppose we observe an 
image using an imperfect apparatus (for example ill-matched glasses), which generates 
some blurring. Such blurring can be described by a so-called point spread function 
specifying the shape generated by a single point of the original image. The observed 
degraded image is consequently given by a convolution of the original image with the 
point spread function. Using this language, we can assign the following names to the 
terms of Eq. (2): 
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W(q',p') original image 

3 — p') 2 ]^j point spread fu 

P(q,p;—\s\) degraded image 



■ exp ( — — -[(g — q') 2 + (p — p') 2 ] ) point spread function 
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The common problem now is the restoration of the original image from the degraded 
one, assuming that the point spread function is known. 



An analytical way to deconvolve Eq. (2) is to apply the Fourier transform, which maps 
a convolution onto a direct product. Dividing both the sides by the Fourier-transformed 
point spread function and evaluating the inverse Fourier transform thus yields an explicit 
expression for the original image. However, this procedure is very sensitive to statistical 
fluctuations and numerical truncation errors, which makes its practical application a 
very delicate matter. These problems have been noted also in the context of quantum 
tomography [9]. 

The numerical instability of the Fourier deconvolution has led to the development 
of techniques dedicated for image restoration. The basic observation is that statistical 
noise docs not allow us to specify precisely the original image that was 'behind' the 
blurred data. In principle, the measured degraded image could could be generated by 
a variety of original images. However, comparing various original images we intuitively 
expect that some of them were more likely to generate the measured data than other 
ones. The maximum-likelihood methodology quantifies this intuition, and selects as an 
estimate the original image which maximizes the likelihood. 

In order to discuss this idea in detail, we shall consider a discretized version of 



where w n is the original image, A vn is the point spread function, and p v is the degraded 
image. Note that this formulation is more general compared to Eq. (2), because it allows 
the point spread function to be of different form for each 'element' of the original image 
indexed with n. The likelihood can be quantified using the function 



which has a rigorous derivation when the degraded image is observed as a histogram 
of events governed by Poissonian statistics. The likelihood function for quantum mea- 
surement has been discussed in Ref. [10], where its maximization has been proposed as 
a method for quantum state estimation. 

In classical imaging, it is natural to assume that w n , as well as A un as a function 
of v for each n, are positive definite distributions with sum equal to one. Under these 
assumptions, it is possible to find the maximum of the likelihood function C via simple 
iteration: 
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Eq. (2): 
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which is the essence of the Richardson algorithm for image restoration [5]. A simple 
heuristic derivation of this algorithm can be found in Rcf . [7] . 

Of course, the maximum-likelihood approach is not a magic wand solving uncon- 
ditionally the problem of image restoration. With increasing blurring, the quality of 
the reconstructed image worsens, and the convergence of the iterative algorithm may 
be very slow. In many cases, however, it offers superior performance compared to the 
Fourier deconvolution technique. 



4 Quantum tomography 

An obvious difficulty with applying the iterative restoration algorithm to quantum to- 
mography is that the object to be reconstructed in the quantum case, i.e. the Wigner 
function, is not positive definite. Nevertheless, there are some other quantum mechan- 
ical reconstruction problems, where positivity constraints appear in a natural way. An 
interesting and nontrivial example is determination of the photon number distribution 
from random phase homodyne statistics [11]. The relation between the phase-averaged 
homodyne statistics and the photon number distribution g n is given by 

J d6h(x;6)=Y / A n (x) gn , (6) 

where A n (x) describe contributions to the homodyne statistics generated by different 
Fock states \n). This formula, after discretization of x, is exactly of the form assumed 
in Eq. (3). Thus we arrive at the following formal analogy: 

g n original image 

A n (x) point spread function 



i r 

— d6 h(x; 9) degraded image 

j7r J -IT 
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which allows us to apply directly the iterative reconstruction algorithm [12]. In this 
procedure, there is one a priori parameter: it is the cut-off of the distribution g n 
specifying the maximum number of photons. 

Reconstruction of the photon distribution may seem to be quite distant from the 
starting point of our considerations, which was deblurring of the Wigner function. How- 
ever, let us recall that the Wigner function can be represented as an alternating sum 
of the photon distribution g n {q,p) corresponding to the phase space displaced state 
D^q,p)gD(q,p): 

1 OO 

W(q,p) = -J2(- 1 ) n Qn(q,p) (7) 

Obviously, we can apply this formula to evaluate the Wigner function at q = p = 0. 
What would be of interest, is the generalization of the maximum-likelihood algorithm to 
determination of an arbitrarily displaced photon distribution g n (q, P )- This would yield 
a numerically stable procedure for reconstructing the Wigner function from homodyne 
statistics, even in the case of the non-unit detection efficiency. 
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Fig. 1. Reconstruction of the Schrodinger cat state |^) oc \2i) — \ — 2i) from Monte Carlo 
simulated homodyne experiment. The homodyne data consisted of 10 5 events generated for 
each of 64 phases spaced uniformly between and 7r. At each point of the grid, the displaced 
photon statistics was obtained from 10 4 iterations, starting from a flat distribution for < 
n < 39. In the simulations, the homodyne variable x has been discretized into 16000 bins over 
the range — 8 < x < 8. 

Surprisingly, this generalization is quite straightforward. The basic observation is 
that the displacement transformation has a simple effect on the homodyne statistics, 
shifting it by y/rj(q cos 9 + p sin 9) for a given local oscillator phase 9. Consequently, we 
have the relation 



which can be readily implemented in the iterative algorithm. Thus, we have arrived 
at the following two-step algorithm for quantum tomography: for a given phase space 
point (q,p), construct the phase-averaged histogram according to the right-hand side 
of Eq. (8), and iteratively reconstruct Qn(q,p)- Then, calculate the value of the Wigner 
function according to Eq. (7) . In Fig. 1 we present Monte Carlo simulated reconstruction 
of the Wigner function for a Schrodinger cat state detected using a homodyne setup 
with the efficiency r\ = 90%. 
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The standard filtered back-projection algorithm used in quantum tomography is 
based on the inverse Radon transform, whose integral kernel is singular. Therefore, a 
regularization scheme is necessary in processing experimental data. This aspect has a 
counterpart in the maximum- likelihood algorithm. In this approach, we have the cut-off 
for the photon distribution which can be regarded as a regularization parameter. Its 
proper choice is an important matter. Setting it too small perturbes the reconstructed 
photon distribution. On the other hand, the larger number of g n s, the slower iterations 
converge. The expected shape of the photon distribution can be quite easily predicted, 
if we roughly know the region of the phase space where the Wigner function is local- 
ized. For this purpose it is useful to recall the semiclassical picture of projections on 
Fock states as rings in the phase space characterized by the radius \/2n. The photon 
distribution is nonzero over the range of n for which the corresponding rings overlap 
with the localization region for the Wigner function. 

Truncation of the photon distribution can be introduced as a regularization scheme 
also in the standard linear reconstruction approach. In such a scheme, the Wigner 
function would be evaluated from a finite part of the photon distribution reconstructed 
using the pattern function technique. However, properties of the reconstructed photon 
distribution make this method very sensitive to statistical noise. This can be straight- 
forwardly seen in the most regular case of r\ = 1. For large n, the error of g n tends 
to a fixed nonzero value [13], and moreover deviations of consecutive g n s are strongly 
anticorrelated. Consequently, the alternating sum defined in Eq. (7) accumulates the 
statistical uncertainty of the photon distribution [14]. 

Let us also note that in principle we could apply the restoration algorithm to homo- 
dyne histograms described by Eq. (1), in order to obtain deblurred quadrature distribu- 
tions (x0\g\xg). In this case statistical fluctuations would play a much more significant 
role. The advantage of using Eq. (8) is that we use the whole available sample of 
experimental data to determine a relatively small number of parameters g n {<l,p)- 
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